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Abstract 

Some biochemical systems show oscillatory behavior. These systems often 
consist of negative feedback loops with repressive transcription regulation. 
Such biochemical systems have distinctive characteristics in comparison with 
ordinary oscillatory chemical systems: i) the numbers of molecules involved 
are small, ii) there are typically only a couple of genes in a cell with a finite 
regulation time scale. Due to the fluctuations caused by these features, the 
system behavior can be quite different from the one obtained by deterministic 
rate equations, because the rate equations ignore molecular fluctuations and 
thus are exact only in the infinite molecular number limit. The molecular 
fluctuations on a free-running circadian system have been studied by Gonze 
et al. (2002) by introducing a scale parameter Vt for the system size. They 
consider, however, only the first effect, assuming that the gene process is fast 
enough for the second effect to be ignored, but this has not been examined 
systematically yet. In this work, we study fluctuation effects due to the finite 
gene regulation time by introducing a new scale parameter r, which we take 
as the unbinding time of a nuclear protein from the gene. We particularly 
focus on the case where the fluctuations due to small molecular numbers 
can be ignored. In simulations on the same system studied by Gonze et 
al., we find the system is unexpectedly sensitive to the fluctuation in the 
transcription regulation; the period of oscillation fluctuates about 30 min 
even when the regulation time scale r is around 30 s, that is even smaller than 
1/1000 of its circadian period. We also demonstrate that the distribution 
width for the oscillation period and the amplitude scales with a/t, and the 
correlation time of the oscillation scales with 1/r in the small r regime. The 
relative fluctuations for the period are about half of that for the amplitude, 
namely, the periodicity is more stable than the amplitude. 
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Introduction 



One of the outstanding features in biological systems is that the systems 
often operate on surprisingly small numbers of active molecules, yet they 
seem to work quite reliably. This is especially intriguing in the case where 
the chemical reaction system involves a gene transcription because there are 
typically only a couple of genes in a cell. 

One example is a circadian system, which shows a rhythmic behavior 
of approximately 24-hour periodicity. It is a universal feature of biological 
systems and know n to be very accurate and rob ust against external and inter- 
nal perturbations ( iDunlapl . Il999l : lYound. 120001). Its biochemical mechanisms 



have been proposed in several systems (IGoldbeterl . Il995l . Il996t iLeloup et al 



19991 ). and most of them are based on a time-delayed negative feedback loop 
of a biochemical reaction network which includes transcription regulations. 
Some of the protein molecules are expected to be very small in number, and 
the number of each gene is typically of the order of one in a cell and does 
not scale with the cell size, thus it is surprising that the circadian system is 
capable of maintainin g its extraordinary regular ity especially in the case of 
a single cell organism ( iBarkai and Leiblerl . 120001 ) . 

The ef fects of mole c ular fl uctuations on the circadian system has been 
studied by Gonze et al. ( 2002bl a) bv Monte Carlo simulations using the Gille- 
spie method f Gillespi 19761 . 1977 ) with the scale parameter for the molec- 
ular numbers. By simulating the system with various values of the scaling 
parameter Q, they demonstrated that the system shows reasonably coherent 
oscillation as long as the system contains more than several tens of mRNA, 
thus concluded that their system are fairly robust against molecular fluctu- 
ations. 

They examined the system rather systematically based up on a standard 
method to study the stochastic nature of chemical reactions (INicolisl . 119721 ) 
by scaling the reaction rates in the way to keep the rate equation unchanged. 
However, there is an ambiguity in the treatment of the gene regulation pro- 
cess because the number of genes should not scale with other protein num- 
bers. They scaled the reaction rates proportional to fl for the gene processes, 
namely, they made the gene regulation times infinitely fast in the large Q 
limit, thus the system dynamics reduce to the one d escribed by t h e corre - 



sponding rate equations without the gene processes (IGonze et al.l . l2002bl ). 
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This could be justified only when the gene processes are so fast that they 
do not cause significant fluctuation on the system behavior. In fact, it is 
not reasonable to assume that the time scale of the gene regulation depends 
upon the scale parameter Q, if you think of it as the cell volume, because 
the time scale with which a regulatory protein binds to the operator site is 
determined by the protein concentration, and the time scale with which the 
protein unbinds is determined by its binding energy. 

In this work, in order to analyze fluctuations from the two distinct ori- 
gins separately, we introduce a new scale parameter r in addition to Q. The 
scale paramter r scales the binding/unbinding time of the gene regulatory 
protein. Thus, these two parameters, r and control the two distinct fluc- 
tuation sources that exist in the biological systems, namely, Q controls the 
fluctuations due to the finite molecular numbers while r controls the fluctu- 
ations due to the finite gene regulation times. We per form the Monte Carl o 
simulations on the same system as the one studied by lGonze et al.l (l2002al ). 
focusing on the latter effect, and demonstrate that significant fluctuation can 
arise from the stochasticity in the gene process alone. We also examine how 
the fluctuation scales with r. 



Model 



The model we study is the simplest version of core model for a circadian 
system that consists of a gene G, mRNA M, cytosolic protein Pc, and nuclear 
protein Pn (FigH]). The biochemical reactions for these elements are given 
by 
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The parameter n is the number of nuclear proteins Pn that bind to suppress 
the gene, i.e. Hill co efficient for the gene activity; we adopt n = 4 for most 
of the calculation as Gonze et al. 
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Figure 1: Simplified core model for a circadian system. 



Now, we introduce the two scaling parameter f2 and r; f2 scales the re- 
action rates so that the numbers of mRNA and the proteins become propor- 
tional to it, and r scales the binding/unbinding time of the nuclear protein 
to the gene operator site. The transition rates for each reaction are listed in 
Table [H where we define the variables G, M, Pc, and P/v to represent the 
numbers of active genes, mRNA, cytosolic proteins, and the nuclear proteins 
in a single cell, respectively. The gene variable G takes either 1 or values, 
depending upon the active state (G) or the inactive state (GPNn), respec- 
tively. Note that we employ Michaelis-Menten enzymatic reactions for the 
degradation processes. The first two reaction rates in Table [1] are for the 
gene regulation and proportional to 1/r, but do not scale with fl because we 
assume only one gene in a celfl. The ratio of the binding and the unbind- 
ing times is determined in the way that the corresponding average behavior 
described by the rate equation remains the same with the original system in 
the small r limit. On the other hand, the gene transcription activity in the 
third reaction is scaled as in order that the numbers of mRNA and the 
proteins should be proportional to Q. 

If we naively write down differential equations for the time evolution. 



^ Gonze et al. ( 2002lJ ) used £7 as a scaling parameter in the places where we use l/r. 
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ignoring the fact that the variables are integers, we would have 

— = vsnc-v^^n^^j^, (8) 
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For ordinary chemical reactions without a gene transcription process, the 
stochastic dynamics should be well described by such equations for the large 
VL case, where the numbers of molecules are large. However, in the present 
system, the number of gene is one and does not scale with fi, thus the 
stochastic nature remains even in the case of the infinite f2 as long as r is 
finite. 

Large Vt limit: 

Supposing the scale parameter as a cell volume, we define the "concen- 
trations" of mRNA and the proteins as, 

[Ml.f, IPcl.f. [P.l.f. (12) 

and write down the rate equations for them as 

%i = *.[M1-..^J^|^-*.[Pc1 + *.|Pk]. (14) 

%i = *i[Po1-*2|Pn1, (15) 

where we have ignored the term of the order oi n/VL in Eq.( |T5l) . 

For ordinary chemical reactions, we expect that the deterministic dynam- 
ics represented by the rate equations would describe the system accurately 
in the large VL limit, because the effect of molecular fluctuation becomes neg- 
ligible. However, for the present case, the system remains stochastic even in 
the large VL limit because the variable G remains stochastic. 
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Small T limit: 

In the case where r is much smaller than a ny other time scales in the 



system, the system reduces to the one studied by lGonze et al.l (j2002al ). This 
can be seen by introducing the time dependent average value of G, denoted 
by Gav(^), i-e. the time average of G over the longer time scale than r 
but shorter than other time scales in the system. Its value is given by the 
condition that the first two reactions in Table [1] are equilibrated, 

Then the system dynamics are given by the stochastic dynamics of reaction 
1 ~ 6 in Table [1] with G replaced by Gav 

If we take the large Q limit on top of this, we obtain the deterministic 
rate equations given by Eqs. f|T3|) ^ f|T5|) with G being replaced by G^vit) of 

Simulations and results 

In order to examine the effect of gene fluctuations, we have performed 
numerical simulations for various values of r and Q. We examine two cases: 
(i) the case where both r and fl are finite, and (ii) the case where r is finite 
but in the large Q limit. In the first case, the fully stochast ic dynam i cs are 



given by Table [U for these we employ the Gillespie method (jGillespid . Il976l . 



19771 ). In the second case, the concentrations [M], [Pc], and [Pn] follow the 
deterministic dynamics while the gene process remains stochastic. In this 
case, we integrate the rate equations ( fT3l) ~ (fT5l) using Runge-Kutta method, 
but at every time step of the length At, the gene variable G is subject to a 
trial for change according to the probability wAt under Poisson process with 
w being the transition rate given in the first two processes in Table 1. 

Figure [2] shows the system behaviors for various values of Q with r = 



O.Olh. The other reaction parameters are the same as those used by lGonze et al 



(j2002al ). The plots in the left column are the time variations of the concen- 
trations of mRNA (solid lines) and the cytosolic protein (dashed lines), and 
those in the right column show the oscillation trajectories projected on the 
M/Q — Pc/^ plane of the phase space. The fluctuation decreases as Q 
increases, but it remains finite even in the infinite Q case because of the 
fluctuation from the gene activity. 
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1= 0.01 h, £2 = 30 



1=0.01 h, £2 = 100 



1= 0.01 h, £2 = 300 



1= 0.01 h, £2 = oo 




Figure 2: Oscillatory behaviors of the concentrations of mRNA and the cytosolic protein 
for various values of fl with r = 0.01 h. The plots in the left column show the time variation 
and those in the right column are the projections of the trajectories in the M/fi — Pc/^ 
plane. We employ the same reaction parameters with those in Gonze et al.; n ~ A, Ve = 0.5 
nM h-i, Ki = 2.0 nM, = 0.3 nM h-\ Km = 0.2 nM, k, = 2.0 h-\ va = 1.5 nM h'^, 
Kd = 0.1 nM, ki ^ 0.2 h-i, k2 = 0.2 h^^. 
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Table 1: Reaction table for a simplified circadian system. 



In order to see the effect of gene stochasticity, we examine the case for 
various values of r in the large Q hmit (Figj3]). The fluctuation decreases on 
decreasing r as in the case of increasing Q. The trajectories are smoother in 
comparison with the previous case because the stochasticity is limited to the 
gene activity. One can see that the fluctuation is evident even in the case 
r = O.Olh, where the ratio r to the period (~ 24h) is as small as 0.5x10"^. 

Figure H] shows the period (i.e. the peak-to-peak interval) distributions 
and the peak value distributions of [Pc] for r = 1.0, 0.1, and 0.01 h. The 
averages and the standard deviations for the distributions are tabulated in 
Table HI Both of the distributions becomes narrower for the smaller value of 
T approximately as a/t , but the standard deviation of the period distribution 
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Figure 3: Oscillatory behaviors of the concentrations of mRNA and the cytosolic protein 
for various values of r with Q ~ 00. The parameters are the same with those in FiglJJ 



is still about a half hour even for the case of r = 0.01 h. It should be noted 
that the ratios of the standard deviation to the average for the peak value 
distributions are about twice as large as those for the period distributions. 
More systematic data are presented in Supplementary material to show the 
a/t scaling and the ratio of the two distribution widths. 

The time correlation function C(t) of the nuclear protein concentration 
[PN(t)] is defined as 

C(A = A[PN(t' + 1)] A[PN(t')] dt', (17) 

where A[PN(t)] represents the deviation from the average, 

A[PN(t)] = [PN(t)]-i^^[PN(t')]rft' (18) 
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(a) Period distribution 



1= 1.0 
0.1 
0.01 




|||Uy_ 



Period (h) 

(b) [Pq] peak distribution 




[Pel (nM) 

Figure 4: Distributions for (a) the period (i.e. peak-to-peak interval) and (b) the peak 
value of the cytosolic protein variation for r = 1.0, 0.1, 0.01 h with il = oo. The averages 
and the standard deviations are tabulated in Table [2l from which one can sec that the 
width of the distribution scales roughly as ^/t. 



r (h) 




1 


0.1 


0.01 




av. 


24.6 


22.3 


22.0 


period (h) 


std. 


3.92 


1.32 


0.45 




std./av. 


0.16 


0.059 


0.020 




av. 


7.82 


6.50 


6.37 


peak value (nM) 


std. 


1.87 


0.79 


0.26 




std./av. 


0.24 


0.12 


0.041 



Table 2: The averages and the standard deviations for the period and the peak value 
distributions for the cytosolic protein concentration shown in Fig|4l The ratios of the 
standard deviation to the average for the peak value distributions are about twice as large 
as those for the period distributions. 
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with T being the time length of the whole simulation. In FigJSl the correlation 
functions are plotted and fitted to the form of damped oscillation 



A cos(a;o^ + ^o) e"*/""^°" (19) 

to estimate the correlation time Tcorr- Figure |6] shows the r dependence of 
the correlation time Tcorr in the logarithmic scale. It shows the scaling 

1 , , 

Tcorr ~ - (20) 

r 

in the small r regime, and the longer correlation time in the n = 4 case than 
in the n = 1 case. One may notice that the correlation time for r = 0.01 h 
is quite long, i.e. Tcorr ~ 2000 h for the n = 4 case, even though the period 
fluctuations are substantial as can be seen in Fig. |4] (See Appendix). 

The scaling of Tcorr given by Eq.( l20|) can be understood as the phase 
diffusion when the standard deviations of the period distributions scales as 
/t as shown in Fig. |H 



Discussion 

We have examined the effects of molecular fluctuations in a biological 
system on a simplified model of a circadian rhythm system, where there are 
two types of fluctuation sources: (i) small numbers of molecules involved 
and (ii) fi nite time sca l e of th e gene regulation. The first effect has been 



studied by lGonze et al.l ( l2002al ). assuming that the gene regulation time scale 
is infinitesimal. In the present work, we focus on the second effect, i.e., in 
the case where the molecular numbers are large enough that the fluctuation 
due to the first effect is neghgible. 

We have developed a method to study this effect systematically by intro- 
ducing a new parameter t to scale the gene regulation times. We set t to be 
the unbinding time of the transcription factor, keeping the ratio of the bind- 
ing to the unbinding rate constant. We performed numerical simulations for 
various values of t without an external entrainment of the 24-hour period. 
As T decreases, the oscillation appears more deterministic; the width of the 
distributions for the oscillation periods and the peak values scales with ^/T 
and the correlation time for the correlations function scales with 1/t. We 
have found that the system is very sensitive to such fluctuation, and demon- 
strated that the oscillation period fluctuates by about 30 min even for very 
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Figure 5: The time correlation functions for the nuclear protein concentration for the 
regulation time r =0.01, 0.1, 1 h with il = oo. The (green) lines shows the fitting curves 
of the form Acos{ujot + Oq) e'*/'^'^°". The fitted values of Tcon- are shown on the plots. The 
other parameters are the same with those in FigI2] 



12 





1 1 




1 1 1 1 

n = 


4 


"1 ' . 

+ 










1 


X 


1000 












T 

^corr 












100 






X 


x^ 


+ 












10 


■ 1 




..1 




..1 , ■ 




0.01 




0.1 


X 


1 



Figure 6: The regulation time t dependences of the correlation time Tcorr for n = 1 and 
4 with = oo in the logarithmic scale. The (green) dashed lines shows the fitted lines 
proportional to t~^. 



small T = 0.01 h 30 s in comparison with its period around 22 li. For the 
present parameter set, the nuclear protein concentration P/v/f2 oscillates in 
the range ~ 5 nM, therefore, the value of r ~ 30 s for the unbinding time 
gives about 1 s for the binding time. These estimates may be tested with 
experimental data. 

The 30 minutes period fluctuation is large for a circadian system. This 
sensitivity to the fluctuation in the gene regulation is an interesting feature 
of the present simplified model. Multiple feedback loops with several phos- 
phorylation steps found in actual bi ological systems m a y be designed in such 



a wa y as to reduce this sensitivity (jUeda et al.l . l200ll : iKusakina and Doddl . 



20121 ). This can be studied by extending the present method. 



The correlation function for the protein oscillation fits to the damped 
sinusoidal function very well, and the estimated correlation time Tcorr scales 
as 1/r in the small r regime. Such decay in the correlation function is 
caused by the phase diffusion due to fluctuations. We estimate Tcorr for the 
Hill coefficient for the gene regulation n =4 and 1, and found that Tcorr 's for 
n = 4 are about 5 times larger than those for n = 1; the fluctuation effect 
is suppressed by the la rger value of t he Hill coefficient by the cooperativity 
effects as in the case of Gonze et al. 

It is also interesting to find that the relative fluctuations for the peak 
values are twice as large as those for the periods, namely, the period is more 
stable than the amplitude. This may be a reason why the correlation time 
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is quite long in spite of apparent fluctuations in the oscillation. 

In the present work, we study only the case where the copy number of the 
gene is 1. However, there are typically a couple of genes in a cell. In the case 
of multiple genes in a cell, the fluctuation in each gene cancels each other, 
therefore overall fluctuation will be reduced. We confirmed by simulation 
that the fluctuation for a two-gene system with r is almost the same as that 
for a single gene system with r/2. This is because the fluctuation cancellation 
by two genes should be comparable with that by one gene that switches twice 
as fast. Simulation data are presented in Supplementary material. 

The fluctuation indicated by our simulations may be compared with pre- 
vious experimental observations. Although circadian clocks are very accu- 
rate as a system, large fluctuatio ns have been o bserve d in the oscillation 



of individual cells of fi broblasts f lNagoshi et al.l . |2004| ) and cyanobacteria 



(IMihalcescu et al.l . |2004| ) when they oscillate independently. For both cases, 
it is reported that the fluctuations are i nuch larger for the amp litude than 
those for the period. In the latter case (IMihalcescu et al.l . 120041 ). the corre- 
lation time is estimated as long as 166±100 days in spite of apparent large 
fluctuations in the amplitude. Such a long correlation time corresponds to our 
case of the gene regulation time scale r = 0.01 h, which gives the r^orr ~ 1950 
h. 

Very little fluctuation is usually observed in circ adian systems; f l uctua - 



tion in the period is typically less than 10 minutes (lAmdaoud et al.l . l2007l ) 



which is even smaller than the fluctuation of 30 minutes that we obtained 
for the case r = 0.01 h. There are some possible mechanisms to suppress 
molecular fluctuations, (i) Cooperativity among cells: The present system 
models a single cell behavior, but for the case of multicellular organisms, 
the cooperativity among cells may exist and that should reduce the fluctu- 
ation in each cell. Actually, variability in each cell is much larger than that 



of a whole system in the case of multicellular organisms fiLiu et all. I1997I : 



Yamaguchi et al.l . l2003l : iHerzog et al.l . l2004l : ICarr and Whitmord . l2005l ) . (ii) 



Multiple feedback loops: Our model is a simplified core model for a circadian 
system and consists of a single negative feedback loop. However, it has been 
known that circadian systems typically consist of multiple feedback loops 



(jZeng et al.l . Il996l : iLeloup and Goldbeterl . Il998l : iGlossop et al.l . Il999t iBlau 



200 ll ). which could be designed in the way to compensate the fluctuations in 
one loop by the other, (iii) Chemical oscillation without gene control: In the 
case of cyanobacteria, it has been proposed that the circ adian system consist s 



of proteins only and does not involve a gene expression (ITomita et al.l . |2005| ) 
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In such a system, the fluctuation discussed in this work does not exist. 

Other than circadian systems, there are some oscillati ons observed i n bio l- 
ogy such as Hesl oscillation during somite segmentation ijHirata et al. . 2002 ). 



p53 o scillation after DNA damage by gamma irradiation fiGeva-Zatorskv et al. 



2000 



2006 ) , or oscillations i n arti f icially constructed s ystems ( lElowitz and Leibler 



Atkinson et all . l2003l : iTigges et al.l . |2009[ ). In these systems, the flue 



tuations are much more profound than circadian systems, and part of the 
fluctuations should come from the gene regulatory processes, for which the 
present analysis is applicable. 

In summary, we have developed a theoretical tool to study the molecular 
fluctuation due to the finite transcription regulation time, and have demon- 
strated that a symplified core model of circadian system is sensitive to such 
fluctuation. Our method can be extended to study a more realistic system 
and can be utilized to clarify biological significance of a detailed design of 
circadian system in terms of stability against the molecular fluctuation. 



Appendix: Phase diffusion and correlation time 



In the appendix, we show that the correlation time Tcorr in the correlation 
function is proportional to 1/r when the period distribution has the width 
proportional to ^r- 

Suppose the correlation function C{t) is written as 



C{t) = A / cos(wot + 0)P{9, t)d9 



(A.l) 



in terms of the average over the phase difference 6 by the distribution func- 
tion P{6,t) at time t. Here, uq is the average angular velocity given by 
Uq = 27r/To in terms of the average period Tq. Now, we assume that the 
phase distribution can be approximated by the Gaussian distribution with 
the standard deviation (Tg{t), 



P{9,t) 



exp 



v/2^^ 
then, Eq. llA.ll) may be estimated as 

C{t) ^ Acos{uot) exp 



2al{t) 



(A.2) 



(A.3) 
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Now, we estimate the phase distribution P(6', t) as follows. The phase 9 
at the time t = uTq may be expressed as the sum of n phases accumulated 
by the time: 

where Tj is the z'th period (i.e. peak-to-peak interval) with Tj = Tq + ATj, 
and we have assumed ATj <C Tq in the last approximation. 
If there is no correlation among ATj, then ag{t) is given by 

alit) ^ (^27r^J (A.5) 

where ay is the standard deviation of the period T. We have replaced n by 
t/To. 

With Eq.f lA.3p . this gives the form of Eq.(IT9]) with 60 = and 

1 T \ ^ 

thus if (Tt oc -\/r ) we obtain Tcorr oc 1/r. 
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Supplementary material to 
Transcription fluctuation eff^ects on biochemical oscillations 

by Ryota Nishino, Takahiro Sakaue, and Hiizu Nakanishi 



In this supplementary material, we present data to demonstrate the a/t- 
scaling of the distribution width and the multiple gene effects more system- 
atically. 

S.l. y^-scaling of the distribution width 

In the text, we show only three sets of data for the distributions of the 
period and the peak value of [Pc] in order to demonstrate the ^ scaling of 
the distribution width. We examined this scaling more systematically. 

In Fig JS.lt the ratios of the standard deviation to the average are plotted 
against r in the logarithmic scale. The dashed lines denote the fitting lines 
with the slope 0.5, which shows the yjr scaling of the distribution width of 
the oscillation parameters in the small r region. 

This can be understood naturally; The on/off frequency of the gene is 
proportional to 1/r, thus the fluctuation in the total on/off-time length scales 
with a/t, from which we expect that the distribution width of the oscillation 
parameters scales with a/t because the response in fluctuation should be 
linear in the small input fluctuation limit. 
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Figure S.l: Scaling of the distribution width by t. The ratio of the standard deviation to 
the average is plotted against r in the logarithmic scale for the peak value of [Pc] and the 
period distributions. The dashed lines are the lines with the slope 0.5 fitted to each data 
set. 



S.2. Multiple gene effects 



In the paper, we only consider the case where the copy number of the gene 
is one. Here, we present the analysis for the case where the copy number of 
the gene is Go- The transition rate for each reaction is given by Table IS.1[ 
where G is the number of active genes. Note that the transition rate for the 
reaction 1 is scaled by I/Gq in order to give the same rate equation as before 
in the limit of r — and Vt — )■ oo. 

In the r — )■ limit, the time dependent average of G denoted by Gx^it) 
is given by 

^^^^^^ ^ ^°l + ([PN]/i^/)"' ^^-^^ 

which gives the same rate equations for the concentrations if the transition 
rate for the reaction 1 is scales by I/Gq. 

To see the copy number effects, we performed Monte Carlo simulations 
for Go = 1 and 2 with some values of r in the f2 — ?■ oo limit. Fig. IS. 21 shows 
the distributions for the period and the peak values of [Pc]- One can see 
the distribution for Gq = 2 with r = 0.1 h agrees with that for Gq = 1 with 
r = 0.05 h quite well for both of the distributions. 

Actually, this can be understood in a simple way; For the copy number 
Go = 2, the fluctuations in the activity of the two genes cancel each other. 
This cancellation should be comparable to the fluctuation cancellation in the 
system with Go = 1 and the half time scale because the gene activity switches 
between on and off twice as fast. 
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Tabic S.l: Reaction tabic for a simplified circadian system in the case where the copy 
number of the gene is Gq . 
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Figure S.2: Period and [Pc] peak value distributions for Gq ~ 1 and 2 with various values 
of T in the 51 — > oo limit. The other parameters are the same with those in the text: n = 4, 
Vs = 0.5 nM h-i, Ki = 2.0 nM, «,„ = 0.3 nM h-\ A',„ = 0.2 nM, = 2.0 h-\ Vd = 1.5 
nM h-i, Kd = 0.1 nM, fci = 0.2 h^^ = 0.2 h"!. 



